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Abstract 

We explore possible applications of the Lattice-Boltzmann Method for the simula¬ 
tion of geophysical flows. This fluid solver, while successful in other fields, is still 
rarely used for geotechnical applications. We show how the standard method can be 
modified to represent free-surface realization of mudflows, debris flows, and in gen¬ 
eral any plastic flow, through the implementation of a Bingham constitutive model. 
The chapter is completed by an example of a full-scale simulation of a plastic fluid 
flowing down an inclined channel and depositing on a flat surface. An application 
is given, where the fluid interacts with a vertical obstacle in the channel. 
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1 Introduction 


Geophysical flows are dangerous natural hazards occurring mostly in moun¬ 
tainous terrain. The most apparent phenomena of this category are debris 
flows, which originate when heavy rainfall mobilizes a large amount of debris 
ng. The resulting mixture comprises water, cohesive sediments, organic mat¬ 
ter, silt, sand and in many cases also stones of different sizes The resulting 
rheological behavior is known to have a wide variability [21], which makes nu¬ 
merical studies an essential tool to support experimental investigations jUj- 
Full-scale simulations of geophysical flows are very scarce, since they require a 
framework that efficiently manages complicated boundary conditions, as well 
as a powerful and flexible fluid solver. Moreover, the non-Newtonian nature 
of the material, and in some cases its multiple phases, pose more challenges. 
Traditional solvers are known to have troubles in tackling this problem, and 
nowadays alternative solution are sought by the community. 

The Lattice Boltzmann Method (LBM) [23| is becoming increasingly popu¬ 
lar and is today considered a valid alternative for categories of flows where 
traditional solvers exhibit disadvantages, like multiphase fluids, flows through 
porous media [IH], irregular geometries ra. and free-surface realizations [T2] . 

After reviewing the most commonly used rheological model for flowing geoma¬ 
terials in Sec.[^ we offer an essential overview of the method in Sec.|^ together 
with a simple but effective formulation for the simulation of geophysical flows. 
In Sec. and examples are given. 


2 Rheology of geophysical flows and simulation 


The rheology of geophysical flow materials is a debated issue in the held, due 
to the extreme variability in natural material parameters and the presence of 
multiple phases, complicating the classification. Most models therefore adopt 
simplified solutions based on single-phase descriptions. This can either be a 
frictional material |2TTT6| or a viscoplastic fluid um The former is used for 
rock and snow avalanches, while the latter is preferred for mudflows and vis¬ 
cous debris flows [2|. For certain categories of geophysical flows, however, a 
single-phase approach is insufficient to capture the physics of the phenom¬ 
ena. Debris flows are a typical example of this, because granular and viscous 
behavior interact, giving rise to unexpected structures and a localization of 
rheological properties [II]- A continuum-continuum coupling for granular and 
fluid phase is possible, but is incapable of capturing the localization of flow 
properties, which is widely recognized to be a key feature of debris flows. 
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A discrete-continuum approach would of course be able to provide a detailed 
description, but development of this sort of coupling has been slowed by its 
demanding computational cost. This is currently challenged, however, by the 
maturity reached by alternative solvers like Smoothed Particle Hydrodynam¬ 
ics, the Material Point Method, or LBM, which are more flexible in managing 
complex boundary conditions than traditional tools like Finite Differences and 
Finite Volumes. In such methods, the granular phase is treated by a separate 
solver and, for this reason, the fluid model can focus on the nature of the 
material. This is the reason behind our choice to adopt LBM with a purely 
viscoplastic rheological law, an approach that can offer: 

• An efficient framework for the simulation of geophysical flows of plastic 
nature, where the complexity of the boundary does not influence the per¬ 
formance. 

• A convenient environment for the coupling with a discrete method. This 
option opens future chances for full realizations of multiphase flows [TH] . 

Regarding the specific rheological law, we adopt the Bingham model, which is 
widely used to describe plastic fluids due to its conceptual simplicity. It reads: 


7 = 0 if fluid does not yield (a < ay), 

< 

a = ay + (ipi'y if fluid flows {a > ay), 


( 1 ) 


where ay and [ipi denote yield stress and plastic viscosity. An analogous way 
to write the law is through an analogy with Newtonian flow. One defines a 
parameter, the apparent viscosity fiapp, which proportionally relates stress and 
rate of strain and is treated as a variable. In the case of a Bingham fluid, fiapp 
takes the form 


^ l^app'i ^ l^app T 


7 


( 2 ) 


where the apparent viscosity Happ (from now on, for simplicity, called viscosity 
/r), diverges when 7 —>■ 0, which will require special care in the solver. We 
are now ready to introduce LBM in the next section, and to incorporate this 
constitutive law in Sec. o 


3 Lattice-Boltzmann Formulation 


LBM has lately emerged as an attractive alternative to traditional fluid solvers, 
mainly due to its high-level performance and the predisposition to paral¬ 
lelization. LBM is also suitable to the solution of problems involving com- 
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plex boundary conditions [T]. It is beyond the scope of this chapter to give a 
complete description of the method. The reader can refer to Refs. |5I23] for 
a comprehensive review. We will focus on the aspects of the formulation that 
need to be modihed in order to successfully reproduce debris flows. 

In LBM, the fluid is described using a distribution function /j and a set of 
discrete velocities Cj. Density p and velocity u of the fluid are computed as 
the first two moments of the distribution function 


I I 

The evolution of /* is governed by the Lattice-Boltzmann equation 


(3) 


/i(x + 5tCi, t + flt) = /i(x, t) + r2i(x, t), (4) 

where Dj is the operator that represents the effects of inter-particle collisions in 
the fluid. A common way to approximate the otherwise complex expression of 
Dj is the Bhatnagar-Gross-Krook operator j3], which relaxes the distribution 
function to a thermodynamic equilibrium /f®. It can be written as 


Dj = 5t 



(5) 


and features a constant, the relaxation time r, which is related to the kinematic 
viscosity of the fluid p as 


r 



P 


(S) 


With this formulation, and with the setting of a coherent lattice |221 , LBM can 
produce realizations of fluid dynamics in analogy to the Navier-Stokes equa¬ 
tions. The method is accurate in the limit of small Mach number, practically 
Umax < O.OlCg with Cg denoting the lattice speed of sound. We will now de¬ 
scribe two additions to the model necessary for the simulation of geophysical 
flows: a non-Newtonian rheology and a free-surface treatment. 


3.1 Non-Newtonian rheology 


The LBM described in the previous section yields, after the Chapman-Enskog 
expansion jl], the Navier-Stokes equation for Newtonian fluids. A simple way 
to upgrade the method to more general formulations is offered by a local 


4 




A 



Fig. 1. Representation of the rheology model employed for plastic fluids. The ap¬ 
proximation of the Bingham model is limited by the maximum and minimum values 
for the relaxation time r imposed by the method. Therefore, also the viscosity /ij 
is limited. 

treatment of the relaxation time r [HEg. Any rheological law that can be 
approximated as 

a = p7, (7) 


with /i = /i(7), is suitable for this approach. The relaxation time can in fact 
be directly related to the viscosity through Equation obtaining ad hoc 
formulations for different rheological laws. The Bingham fluid, for example, 
can be written as 


5t 1 ( ay\ 

a = ay + fipa ^ T = — + — I /ip/ + y I • 


( 8 ) 


This type of formulation requires the computation of the shear rate tensor, 
which can be done easily in LBM directly form the distribution functions 

iab = {fi - /D , (9) 

and the magnitude can be extracted as 


^=,/2i;E7a»7.». (10) 

y a b 


The limitation of this approach lies in the range of values given to the re¬ 
laxation time T by Equation Accuracy in LBM is guaranteed as long as 
Tmin < T < Tmax- Reasonable values for these limits are Tmin = 0.501 and 
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Tmax = 1-0. Therefore, also the viscosity jj,, which is linearly linked to the 
relaxation time, is subjected to the same restrictions: fimin < < fJ^max- The 

following considerations are thus necessary: 

• The fluid that reaches the maximum allowed value of /i is considered to be 
in a plastic state. However, with the proposed scheme, the fluid never stops 
its motion, but rather flows at a much slower rate. The ratio between Umax 
and iimin determines the effectiveness of this approach. With the proposed 
limit values for r, fimax = 500fimin- 

• The best approximation of a Bingham fluid is obtained when fimin < 
because the lower limitation on /i has no effect. However, an eventual tran¬ 
sition to turbulent regime can happen when simulating diluted flows, and 
therefore the value of must be raised to avoid instabilities. In case 
fJ'min > IJ'ph the approximation of the Bingham constitutive model becomes 
less accurate. 


3.2 Implementation of the free-surface technique 


In order to simulate geophysical flows on realistic geometries, we need to in¬ 
clude the boundary conditions given by the channel bed and the interface of 
the flow with air. While the former can be implemented as a standard no- 
slip boundary condition, as in Ref. [H], the latter is a less common practice 
in LBM. The free-surface is represented through a classification of the lat¬ 
tice nodes in three categories: liquid, interface and gas nodes. The governing 
parameter is the liquid fraction A: 


X = p if the node is liquid, 

< 0 < A < p if the node is interface, (11) 

A = 0 if the node is gas. 

The liquid fraction of a node evolves according to the streaming of the distri¬ 

bution function given by Equation as 


X(t St) — X{t) H X/ “ /out) (12) 

P 

where fin and font represent the distribution function streaming respectively 
in and out of the node, and a is a parameter that depends on whether the 
distributions are exchanged with a fluid node or another interface node. This 
method conserves mass exactly, and ensures a smooth evolution of the surface. 
Further details are found in Ref. [12]. 
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4 Simulation of mudflow 



Fig. 2. Geometry of the simulation. The fluid mass lies at the top of a long cylindrical 
chute. The deposition area at the bottom is flat and features a vertical obstacle. 
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Fig. 3. Evolution of the maximum velocity of the fluid and of the flow plasticization 
level, computed as ratio between cells that have reached the maximum viscosity and 
the total number of cells. 

The full-scale simulation of a plastic geophysical flow is shown in this sec¬ 
tion. Mimicking the real geometry of a small valley, the simulation features a 
cylindrical channel inclined at 5° with respect to the horizontal, and a flat de¬ 
position area at its bottom (Fig. [^. The total volume of the flowing material 
is 500 m^ and is fixed, i.e. neither entrainment nor deposition are modeled. 
While very big events can be of the order of 10^ m^, the size of the most fre¬ 
quent type of geophysical flows lies in the range of 10^ m^, which is big enough 
to endanger humans and infrastructures. Therefore our simulation proposes 
a realistic scenario, even though not a particularly dangerous one. The fluid 
has density p = 2000 kg/m^ and follows a Bingham-like rheological law, like 
the one proposed in Sec. Yield stress and kinematic plastic viscosity are 
respectively ay = 150 Pa and fipi = 10 textrmrri^/s, relating the simulated 
system to a very dense mudflow or to a debris flow whose granular phase has 
been homogenized into the fluid, therefore increasing the bulk viscosity pnin| . 
Fig^a shows how the fluid free surface evolves in time. The fluid is quickly 
sheared by the effect of gravity and moves until an equilibrium is reached in 
the deposition area, where the viscosity increases. This technique can be used 
to estimate the deposition area of the material after an event, and to support 
the design of hazard maps on real terrain. Fig. shows the evolution of the 
maximum velocity in the fluid and of the plasticization level of the material, 
which are the useful parameters to determine the status of the flow. 
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Fig. 4. Evolution of the geometry of the flow. Intensities show the viscosity at the 
surface, therefore indicating the rate of shearing of the fluid: low (dark) or high 
(light). 

5 Obstacle Interaction 


To show the possibilities to use LBM to design protection structures, we repeat 
the simulation of the previous section, this time featuring an obstacle. LBM 
can in fact be used to calculate the hydrodynamic interactions on solid objects, 
computing all momentum transfers between the distribution function and the 
solid boundaries. The procedure, which is found in Ref. na, does not change 
significantly the overall efficiency of the method. We add a retaining wall, fixed 
at the bottom of the channel and of size H x L x S = 3.0 m x 2.0 m x 0.25 m, 
as in Fig. The shape of the free surface after the impact is shown in Fig. 
with insight into the longitudinal cross section of the flow. 








































Fig. 5. Image of the flow splashing on the retaining wall, at t = 14.0 s. The color 
contour shows the velocity at the free surface (a) and in the longitudinal section 
(b). 
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Fig. 6. Force exerted on the obstacle by the flow. The estimation is obtained with 
the hydrodynamic formula in Eq. 


The force on the wall can be estimated with a hydrodynamic formula [8j as 


Fyoall fronts 


(13) 


where A is the area of the obstacle impacted by the flow. The value of the 
coeflicient k is given by the comparison with experiments and varies, according 
to different authors, from 2 to 5. In the simulation, when the flow hits the 
wall, the depth is 0.5m and the front speed is vfront — 5m/s, which leads to 
an estimated force of Fy^jaii = 300 -i- 750 kN. Fig. shows the hydrodynamic 
force as calculated by the solver, highlighting the importance of the dynamic 
load due to the initial impact. The maximum values match the prediction of 
the hydrodynamic formula. 
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6 Outlook 


In this chapter we showed how a model based on LBM can be used to simulate 
geophysical flows and provide a new tool for the rational design of mitigation 
and protection structures. The model inherits the advantages of the local 
solution mechanism of LBM, and extends the standard solver with the addition 
of a Bingham fluid formulation and of the free-snrface technique. The resulting 
framework can be used to simulate homogeneous plastic flows, and provides 
an optimal environment for the coupling with discrete method, thus opening 
future chances for the full simulation of multiphase geophysical flows. 
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